Calibration Update

library(gotmtools)
## Loading required package: rLakeAnalyzer
## Warning: replacing previous import 'stats::filter' by 'dplyr::filter' when
## loading 'gotmtools'
library(LakeEnsemblR)
## 
## 
##   _          _        _____                          _     _ ____  
##  | |    __ _| | _____| ____|_ __  ___  ___ _ __ ___ | |__ | |  _ \ 
##  | |   / _` | |/ / _ |  _| | '_ \/ __|/ _ | '_ ` _ \| '_ \| | |_) |
##  | |__| (_| |   |  __| |___| | | \__ |  __| | | | | | |_) | |  _ < 
##  |_____\__,_|_|\_\___|_____|_| |_|___/\___|_| |_| |_|_.__/|_|_| \_\
##                                                                    
## 
##               https://github.com/aemon-j/LakeEnsemblR
## LakeEnsemblR version 1.0.8 (2021-08-12) is loaded
## WARNING: Your current version of LakeEnsemblR (v1.0.8) is ahead of the master branch version (1.0.5)
## Development version may have an unexpected behaviour
## 
## Attaching package: 'LakeEnsemblR'
## The following object is masked from 'package:gotmtools':
## 
##     analyse_strat
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(rLakeAnalyzer)
library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:dplyr':
## 
##     rename
library(RColorBrewer)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:reshape':
## 
##     stamp
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/")

ncdf <- "output/ensemble_output_all_models_14Aug21.nc"

model <- c("FLake", "GLM", "GOTM", "Simstrat", "FLake")

spin_up <- 190



plot_heatmap(ncdf, model = model) +
  scale_colour_gradientn(limits = c(0, 32),
                         colours = rev(RColorBrewer::brewer.pal(11, "Spectral"))) + theme_classic()
## Warning in check_models(model): Model: "FLake" is redundant in the input to
## argument "model"
## Extracted temp from output/ensemble_output_all_models_14Aug21.nc
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

plot_ensemble(ncdf, model = model, var = "ice_height")
## Warning in check_models(model): Model: "FLake" is redundant in the input to
## argument "model"

fit <- calc_fit(ncdf, model = model, spin_up = spin_up)
## Warning in check_models(model): Model: "FLake" is redundant in the input to
## argument "model"
fit
## $FLake
##       rmse       nse         r       bias      mae      nmae
## 1 2.337636 0.9050119 0.9531674 -0.1983942 1.496257 0.3488843
## 
## $GLM
##       rmse       nse         r      bias      mae      nmae
## 1 1.951257 0.9314362 0.9716791 0.5829078 1.511608 0.8308405
## 
## $GOTM
##      rmse       nse         r       bias     mae      nmae
## 1 1.85635 0.9379438 0.9700531 -0.3822369 1.26915 0.3418291
## 
## $Simstrat
##       rmse       nse         r      bias      mae      nmae
## 1 1.502343 0.9593553 0.9814047 0.1883044 1.157387 0.2200856
## 
## $ensemble_mean
##       rmse       nse         r       bias       mae      nmae
## 1 1.283356 0.9703407 0.9856149 0.06428875 0.9322176 0.3330117
# out <- analyze_ncdf(ncdf, model, spin_up = 190)
# out$stats 

## Plot residuals
plist <- plot_resid(ncdf = "output/ensemble_output_all_models_14Aug21.nc", var = "temp")
## Extracted temp from output/ensemble_output_all_models_14Aug21.nc
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
ggarrange(plotlist = plist)

Schmidt Stability

  • Not much to say about Schmidt Stability, observations seem to track well with the mean at the very least, RMSE calculations coming soon for the individual models.
library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

ncdf <- "~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc"


######################## Calculating Schmidt Stability ################################


out <- load_var(ncdf = ncdf, var = "temp")
## Extracted temp from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
# out <- as.data.frame(out[[1]])
bathy <- read.csv('~/Dropbox/sunapee_LER_projections/LER_inputs/sunapee_hypso.csv')
colnames(bathy) <- c("depths", "areas")
ts.sch <- lapply(out, function(x) {
  ts.schmidt.stability(x, bathy = bathy, na.rm = TRUE)
})
## Reshape to data.frame
df <- melt(ts.sch, id.vars = 1)
colnames(df)[4] <- "model"
df$yday <- yday(df$datetime)
df$year <- year(df$datetime)

df <- df %>% 
  dplyr :: group_by(yday, year) %>% 
  dplyr :: mutate(mean = mean(value, na.rm = TRUE)) %>% 
  dplyr :: mutate(sd = sd(value, na.rm = TRUE))



## plot results
ggplot(df, aes(yday, value, colour = model)) +
  # geom_line(data=df, aes(y=mean, x=yday), color = "black", size = 1) +
  geom_line() +
  # geom_ribbon(data = df, aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6, 
  #             linetype = 0.1, 
  #             color = "grey") +
  facet_wrap(~year) +
  labs(y = "Schmidt stability (J/m2)") +
  theme_classic() + ylim(-50, 1000)
## Warning: Removed 82 row(s) containing missing values (geom_path).
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggplot(df, aes(yday, mean)) + 
  geom_line() + 
  geom_ribbon(data = df, aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6,
              linetype = 0.1,
              color = "grey") + 
  facet_wrap(~year) + 
  labs(y = "Schmidt stability (J/m2)") +
  theme_classic() + ylim(-50, 1000) + 
  geom_line(data = subset(df, model == "Obs"), aes(yday, value, col = "Obs"))
## Warning: Removed 80 row(s) containing missing values (geom_path).

Thermocline Depth

  • Similar update on thermocline depth. RMSE calculations coming soon.
  • GLM appears to be the poorest performer most years based on the first figure.
  • The extra years of observations give me more confidence, as 2009 and 2010 were a bit off. The thermocline seems to fit quite well from 2011 and on.

All other metrics

  • Observations are added for this group of metrics
  • There seems to be an issue with 2009 calculations, as shown by the printed dataframe.
  • The only metric to be compared with ice will have to be added in - coming soon.
  • There is very limited stratification data calculated for the observed data.
  • Seems to very high uncertainty for bottom temperature. Reason to remove it from desired metrics?
library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

ncdf <- "~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc"

out <- load_var(ncdf = ncdf, var = "temp")
## Extracted temp from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
temp <- load_var(ncdf, "temp")
## Extracted temp from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
ice <- load_var(ncdf, "ice_height")
## Extracted ice_height from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
##            [,1]                  
## short_name "ice_height"          
## units      "m"                   
## dimensions "lon lat model member"
out <- lapply(1:length(temp), function(x) {
  # x = 1 # for debugging
  mlt <- reshape::melt(temp[[x]], id.vars = 1)
  mlt[, 2] <- as.numeric(gsub("wtr_", "", mlt[, 2]))
  if(names(out)[x] == "Obs") {
    analyse_strat(data = mlt)
  }
  analyse_strat(data = mlt, H_ice = ice[[x]][, 2])
})
## Warning: Using 11.5 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 0.5 as the surface.
## Warning: Using 10.5 as the bottom.
## Warning: Using 0.5 as the surface.
## Warning: Using 10.5 as the bottom.
names(out) <- names(temp)


df <- melt(out[1:6], id.vars = 1)
colnames(df)[4] <- "model"


df <- df %>% 
  dplyr :: group_by(year, variable) %>% 
  dplyr :: mutate(mean = mean(value, na.rm = TRUE)) %>% 
  dplyr :: mutate(sd = sd(value, na.rm = TRUE))



ggplot(df, aes(x = year, y = value, colour = model)) + geom_line() + 
  facet_wrap(~variable, scales = "free_y")
## Warning: Removed 17 row(s) containing missing values (geom_path).

ggplot(df, aes(x = year, y = mean)) + geom_line() + 
  facet_wrap(~variable, scales = "free_y") + 
  geom_ribbon(data = df, aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6,
                                                         linetype = 0.1,
                                                         color = "grey") + 
  geom_line(data = subset(df, model == "Obs"), aes(year, value, col = "Obs"))
## Warning: Removed 5 row(s) containing missing values (geom_path).

checkthis <- filter(df, year == 2009, model == "Obs")
checkthis
## # A tibble: 24 x 6
## # Groups:   year, variable [24]
##     year variable      value model     mean      sd
##    <dbl> <fct>         <dbl> <chr>    <dbl>   <dbl>
##  1  2009 TsMax         0.140 Obs    20.8     10.2  
##  2  2009 TsMaxDay      5     Obs   192.      91.9  
##  3  2009 TsMin        -0.270 Obs    -0.0875   0.136
##  4  2009 TsMinDay      0     Obs    74.8    139.   
##  5  2009 TbMax         0.970 Obs     9.46     6.46 
##  6  2009 TbMaxDay     11     Obs   186.     140.   
##  7  2009 TbMin         0.400 Obs     1.27     1.30 
##  8  2009 TbMinDay      0     Obs    34.5     44.5  
##  9  2009 MaxStratDur  NA     Obs   176.      30.8  
## 10  2009 MeanStratDur NA     Obs   129.      65.1  
## # … with 14 more rows

```